/*
	knmsuvival - (c) 2011 Rafael Luiz Klaser

	This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

// knmsurvival.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"

//using namespace std;

#define MAX_ROWS	256
#define MAX_COLS	256

#define ASTAR_COL	0
#define DSTAR_COL	1
#define BSTAR_COL	2
#define W_COL		3
#define A_COL		4
#define SUMA_COL	5

#define NUM_RESULTS	6

char* result_labels[NUM_RESULTS] = {
	"A*",
	"D*",
	"B*",
	"W",
	"a",
	"[a"
};

typedef
struct s_context {
	int timestep;
	int rows;
	int cols;
	int input_matrix[MAX_ROWS][MAX_COLS];
	int astar_matrix[MAX_ROWS][MAX_COLS];
	int astar_matrix_sum[MAX_COLS];
	int dstar_matrix[MAX_ROWS][MAX_COLS];
	int dstar_matrix_sum[MAX_COLS];
	double result_matrix[MAX_COLS][NUM_RESULTS];
	double phi;	//instantaneous mortality rate

}
context;

void load_data(FILE* f, context& ctx);
void calc_astar_matrix(context& ctx);
void calc_dstar_matrix(context& ctx);
void calc_astar(context& ctx);
void calc_bstar(context& ctx);
void calc_dstar(context& ctx);
void calc_phi(context& ctx);
void calc_w(context& ctx);
void calc_a(context& ctx);
void calc_sum_a(context& ctx);
void print_result(context& ctx);
void dump_data(context& ctx);

int main(int argc, char* argv[])
{
	context ctx;
	FILE *fi = NULL;
	int ret = -1;
	

	if(argc<3) {
		std::cout << "Kiritani-Nakasuji-Manly survival estimating method" << std::endl;
		std::cout << "--------------------------------------------------" << std::endl;
		std::cout << "usage: knmsurvival timestep rows cols datafile" << std::endl;
		std::cout << "where," << std::endl;
		std::cout << "\ttimestep is the constant time difference between each data row" << std::endl;
		std::cout << "\trows is the number of data rows" << std::endl;
		std::cout << "\tcols is the number of data columns" << std::endl;		
		std::cout << "\tdatafile is the file containing the data" << std::endl;
		std::cout << "note: first columns is the time column." <<std::endl;
	}
	else {

		memset(&ctx, 0, sizeof(ctx));

		ctx.timestep = atoi(argv[1]);
		ctx.rows = atoi(argv[2]);
		ctx.cols = atoi(argv[3]);

		if(ctx.cols>MAX_COLS || ctx.rows>MAX_ROWS) {
			std::cout << "Oops... we only support " << MAX_ROWS << " rows and " << MAX_COLS << " columns." << std::endl;
			std::cout << "aborting..." << std::endl;
		}
		else {

			fopen_s(&fi, argv[4], "r");
			
			load_data(fi, ctx);
			dump_data(ctx);
			calc_astar_matrix(ctx);
			calc_dstar_matrix(ctx);
			calc_astar(ctx);
			calc_dstar(ctx);
			calc_bstar(ctx);
			calc_phi(ctx);
			calc_w(ctx);
			calc_a(ctx);
			calc_sum_a(ctx);

			print_result(ctx);

			fclose(fi);
		}

		ret = 0;
	}
	
	return ret;
}

void load_data(FILE* f, context& ctx) 
{
	for(int i=0; i<ctx.rows; i++) {
		for(int j=0; j<ctx.cols; j++) {
			fscanf_s(f, "%d\t", &ctx.input_matrix[i][j]);
		}
	}
}

void print_result(context& ctx) 
{
	for(int j=0; j<NUM_RESULTS; j++) {
		std::cout << result_labels[j] << "\t"; 
	}
	std::cout << std::endl;

	for(int i=0; i<ctx.cols-1; i++) {
		for(int j=0; j<NUM_RESULTS; j++) {
			//std::cout << ctx.result_matrix[i][j] << "\t"; 
			printf("%lf\t", ctx.result_matrix[i][j]);
		}
		std::cout << std::endl;
	}
}

void dump_data(context& ctx) 
{
	for(int i=0; i<ctx.rows; i++) {
		for(int j=0; j<ctx.cols; j++) {
			printf("%d\t", ctx.input_matrix[i][j]);
		}
		std::cout << std::endl;
	}
}

/**
 * Cumulative Stage Distribution Matrix
 */
void calc_astar_matrix(context& ctx) 
{
	memcpy(&ctx.astar_matrix, &ctx.input_matrix, sizeof(ctx.astar_matrix));
	for(int i=0; i<ctx.rows; i++) {
		for(int j=ctx.cols-1; j>0; j--) {
			ctx.astar_matrix[i][j] += ctx.astar_matrix[i][j+1];
		}
	}
}

/**
 * Cumulative Stage Frequency Matrix
 */
void calc_dstar_matrix(context& ctx)
{
	memcpy(&ctx.dstar_matrix, &ctx.astar_matrix, sizeof(ctx.dstar_matrix));
	for(int i=0; i<ctx.rows; i++) {
		for(int j=1; j<ctx.cols; j++) {
			ctx.dstar_matrix[i][j] = ctx.dstar_matrix[i][j] * ctx.dstar_matrix[i][0];
		}
	}
}

/**
 * Total Cumulative Distribution per Stage
 */
void calc_astar(context& ctx)
{
	for(int i=0; i<ctx.rows; i++) {
		for(int j=1; j<ctx.cols; j++) {
			ctx.astar_matrix_sum[j] += ctx.astar_matrix[i][j];
		}
	}

	for(int j=1; j<ctx.cols; j++) {
		ctx.result_matrix[j-1][ASTAR_COL] = ctx.astar_matrix_sum[j] * ctx.timestep;
	}
}

/**
 * Total Cumulative Frequency per Stage
 */
void calc_dstar(context& ctx)
{
	for(int i=0; i<ctx.rows; i++) {
		for(int j=1; j<ctx.cols; j++) {
			ctx.dstar_matrix_sum[j] += ctx.dstar_matrix[i][j];
		}
	}

	for(int j=1; j<ctx.cols; j++) {
		ctx.result_matrix[j-1][DSTAR_COL] = ctx.dstar_matrix_sum[j] * ctx.timestep;
	}
}

/**
 * Concentration Time per Stage
 */
void calc_bstar(context& ctx)
{
	for(int j=0; j<ctx.cols-1; j++) {
		ctx.result_matrix[j][BSTAR_COL] = ctx.result_matrix[j][DSTAR_COL] / ctx.result_matrix[j][ASTAR_COL];
	}
}

/**
 * Instantaneous Mortality Rate
 */
void calc_phi(context& ctx)
{
	double a;
	double b;

	a = ctx.result_matrix[ctx.cols-2][ASTAR_COL] / ctx.result_matrix[0][ASTAR_COL];
	b = ctx.result_matrix[ctx.cols-2][BSTAR_COL] - ctx.result_matrix[0][BSTAR_COL];
	
	ctx.phi = -log(a)/b;
}

/**
 * Stage-Specific Mortality
 */
void calc_w(context& ctx)
{
	for(int j=0; j<ctx.cols-2; j++) {
		ctx.result_matrix[j][W_COL] = ctx.result_matrix[j+1][ASTAR_COL] / ctx.result_matrix[j][ASTAR_COL];
	}
}

/**
 * Stage Duration
 */
void calc_a(context& ctx)
{
	for(int j=0; j<ctx.cols-2; j++) {
		ctx.result_matrix[j][A_COL] = -log(ctx.result_matrix[j][W_COL]) / ctx.phi;
	}
}

/**
 * Accumulated Stage Duration
 */
void calc_sum_a(context& ctx)
{
	ctx.result_matrix[0][SUMA_COL] = ctx.result_matrix[0][A_COL];
	for(int j=1; j<ctx.cols-2; j++) {
		ctx.result_matrix[j][SUMA_COL] = ctx.result_matrix[j-1][SUMA_COL] + ctx.result_matrix[j][A_COL];
	}
}
